clear all
clc
close all

load AALabData
Row=row;
Column=column;

load AANewData
Row=[Row;row'];
Column=[Column;column'];

load AAMTData
Row=[Row;row];
Column=[Column;column];

ROW=zeros(486,9);
for i=1:486
    ROW(i,:)=reshape(Row{i},[9,1])';
end
[~,I]=unique(ROW,'rows');

Row=Row(I,:);
Column=Column(I,:);

n=length(Row);

testtau=0:0.001:0.1;
u=length(testtau);

Q=100;
fakedistr=cell(1,Q);
ErrorUniform=zeros(1,Q);
ErrorPCHM=zeros(1,Q);
complete=9999*ones(1,Q);
fittau=zeros(1,Q);

load FM

for q=1:Q
    q
    
    % Generate fake distributions
    fakedistr{q}=genfake(Row,Dominant,Dominated);
    %fakedistr{q}=genfakeNoRestr(Row);
    Freq=fakedistr{q};
        
        % Predict uniform
            
        err=zeros(1,n);
        for i=1:n
            err(i)=Freq(i,1)*log((Freq(i,1)+eps)/(1/3))+Freq(i,2)*log((Freq(i,2)+eps)/(1/3))+Freq(i,3)*log((Freq(i,3)+eps)/(1/3));
        end
        ErrorUniform(q)=mean(err);
        
        % Run PCHM     
        
        levR{1}=1/3*ones(1,3);
        levC{1}=1/3*ones(1,3);
        
        tauerr=zeros(1,u);
        for k=1:u
            tauerr(k)=CHfixedtau(testtau(k),Row,Column,Freq);
        end
        fittau(q)=find(tauerr==min(tauerr));
        ErrorPCHM(q)=min(tauerr);
        
end

'Restrictiveness'
r = mean(ErrorPCHM)/mean(ErrorUniform)

'Standard Error'
sampleVarPCHM = mean((ErrorPCHM-mean(ErrorPCHM)).^2);
sampleVarUniform = mean((ErrorUniform-mean(ErrorUniform)).^2);
covar = mean((ErrorPCHM-mean(ErrorPCHM)).*(ErrorUniform-mean(ErrorUniform)));

SE = sqrt((sampleVarPCHM - 2 * r * covar + r^2 * sampleVarUniform) / mean(ErrorUniform)^2) * (1/sqrt(Q))

